Sequential deterministic optimization based control system and method

ABSTRACT

The embodiments described herein include one embodiment that a control method including executing an infeasible search algorithm during a first portion of a predetermined sample period to search for a feasible control trajectory of a plurality of variables of a controlled process, executing a feasible search algorithm during a second portion of the predetermined sample period to determine the feasible control trajectory if the infeasible search algorithm does not determine a feasible control trajectory, and controlling the controlled process by application of the feasible control trajectory.

CROSS-REFERENCE TO RELATED APPLICATION

Under 35 U.S.C. § 120, this application is a continuation of U.S. patentapplication Ser. No. 13/836,701, filed Mar. 15, 2013, which isincorporated by reference herein in its entirety for all purposes.

BACKGROUND

The invention relates generally to control systems and more particularlyto deterministic optimization based control of systems.

Generally, control system, such as an industrial plant or a powergeneration system, may be dynamic and include various constraints. Forexample, the constraints on the control system may be the result ofactuator limits, operational constraints, economical restrictions,and/or safety restrictions. Accordingly, control of such a multivariableconstrained dynamic system may be complex. Techniques such as coupledmulti-loop proportional-integral-derivative (PID) controllers may not bebest suited for handling the control of such complex control systems. Onthe other hand, one process control technique capable of handling themultivariable constraints is optimization based control (OBC).Specifically, OBC may improve the performance of the control system byenabling the system to operate closer to the various constraints (i.e.,via dynamic optimization).

However, OBC may be computationally demanding because the dynamicoptimization calculation may involve solving a constrained optimizationproblem such as quadratic programming (QP) problems at each samplingtime. Utilizing a general solver may take seconds or even minutes. Inaddition, it may be difficult to predict the time it takes for theoptimizer to solve the constrained optimization problems. Accordingly,to utilize OBC to control systems with faster dynamics, it may often bebeneficial to enable deterministic OBC to provide a feasible controlaction within a predetermined control time.

SUMMARY

Certain embodiments commensurate in scope with the originally claimedinvention are summarized below. These embodiments are not intended tolimit the scope of the claimed invention, but rather these embodimentsare intended only to provide a brief summary of possible forms of theinvention. Indeed, the invention may encompass a variety of forms thatmay be similar to or different from the embodiments set forth below.

A first embodiment provides a control method including executing aninfeasible search algorithm during a first portion of a predeterminedsample period to search for a feasible control trajectory of a pluralityof variables of a controlled process, executing a feasible searchalgorithm during a second portion of the predetermined sample period todetermine the feasible control trajectory if the infeasible searchalgorithm does not determine a feasible control trajectory, andcontrolling the controlled process by application of the feasiblecontrol trajectory.

A second embodiment provides a control method, including executing aninfeasible search algorithm during a first portion of a predeterminedsample period to search for a feasible control trajectory of a pluralityof variables of a controlled process, if the infeasible search algorithmdoes not determine a feasible control trajectory during the firstportion of the sample period, making a current control trajectory of theinfeasible search algorithm feasible and stabilizing at the end of thefirst portion of the sample period, executing a feasible searchalgorithm during a second portion of the sample period to determine thefeasible control trajectory based upon the stabilized current controltrajectory, and controlling the controlled process by application of thefeasible control trajectory.

A third embodiment provides a control system, including memory circuitryfor storing executable code, and processing circuitry for executing thecode. The code defining steps that, when executed executes an infeasiblesearch algorithm during a first portion of a predetermined sample periodto search for a feasible control trajectory of a plurality of variablesof a controlled process, if the infeasible search algorithm does notdetermine a feasible control trajectory during the first portion of thesample period, executes a feasible search algorithm during a secondportion of the sample period to determine the feasible controltrajectory, and controls the controlled process by application of thefeasible control trajectory.

BRIEF DESCRIPTION OF THE DRAWINGS

These and other features, aspects, and advantages of the presentinvention will become better understood when the following detaileddescription is read with reference to the accompanying drawings in whichlike characters represent like parts throughout the drawings, wherein:

FIG. 1 depicts a block diagram of an embodiment of a control systemutilizing deterministic optimization-based control (OBC);

FIG. 2 depicts a block diagram of an embodiment of the deterministicoptimization-based control from FIG. 1;

FIG. 3 depicts a block diagram of an embodiment of the processingcircuitry from FIG. 2;

FIG. 4 depicts an embodiment of a deterministic OBC modulecommunicatively coupled to other modules within the control system.

FIG. 5A depicts an embodiment of a feasible search method for solving aquadratic programming problem that has found an optimum solution;

FIG. 5B depicts an embodiment of a feasible search method for solvingthe quadratic programming problem that has not found the optimumsolution;

FIG. 5C depicts an embodiment of an infeasible search method for solvingthe quadratic programming problem that has found the optimum solution;

FIG. 6 is a flow chart depicting an embodiment of a process controlutilizing sequential solvers;

FIG. 7 depicts a block diagram of an embodiment of a stabilizationprocess used in the optimization-based control;

FIG. 8A depicts a first feasible solution for the quadratic programmingproblem;

FIG. 8B depicts the first feasible solution advanced one step comparedto a second feasible solution for the quadratic programming problem;

FIG. 9 depicts an embodiment of dynamic optimization-based control (OBC)utilizing a steady state block, a dynamic target generator block, and adynamic optimization block;

FIG. 10 depicts an embodiment of dynamic optimization-based control(OBC) utilizing a steady state block and a dynamic target generatorblock;

FIG. 11 depicts a block diagram of an embodiment of the optimizationbased control with a sequential optimization search running in parallelwith another optimization search;

FIG. 12 A depicts an embodiment of complex constraints; and

FIG. 12B depicts and embodiments simple bound constraints acting as afunnel for the complex constraints.

DETAILED DESCRIPTION

One or more specific embodiments of the present invention will bedescribed below. In an effort to provide a concise description of theseembodiments, all features of an actual implementation may not bedescribed in the specification. It should be appreciated that in thedevelopment of any such actual implementation, as in any engineering ordesign project, numerous implementation-specific decisions must be madeto achieve the developers' specific goals, such as compliance withsystem-related and business-related constraints, which may vary from oneimplementation to another. Moreover, it should be appreciated that sucha development effort might be complex and time consuming, but wouldnevertheless be a routine undertaking of design, fabrication, andmanufacture for those of ordinary skill having the benefit of thisdisclosure.

When introducing elements of various embodiments of the presentinvention, the articles “a,” “an,” “the,” and “said” are intended tomean that there are one or more of the elements. The terms “comprising,”“including,” and “having” are intended to be inclusive and mean thatthere may be additional elements other than the listed elements.

The present disclosure is generally directed toward systems and methodsfor deterministic optimization-based control (OBC) of a control system,such as an industrial plant, a power generation system, or the like.Generally, control systems may utilize process control techniques tocontrol the system. For example, some control systems utilizeproportional-integral-derivative (PID) controllers coupled in amulti-loop configuration. Multi-loop PID controllers may offer a fastreal-time control of a control system. In addition, PID controllers mayrun on embedded systems with less computational power. However, whencontrol system has complex dynamics and/or its operation is constrained,the complexity of the process control may greatly increase andmulti-loop PID controllers may not provide adequate control. Forexample, the control system may include processes with large dead timesor non-minimum phase dynamics.

One process control technique for control of dynamic multivariablesystems is optimization-based control (OBC), which can offer bettercontrol (e.g. reduces process variations to enable operation closer toconstraints at more profitable operating points). Specifically, OBC usesa process model to predict future process output trajectories based onprocess input trajectories. In other words, OBC computes trajectories ofmanipulated variables to optimize the objective function (i.e., minimizecosts). As used herein, the cost includes the determination of how wellthe output trajectories match the desired setpoints. It should beappreciated that in linear control systems, the cost may be captured asa quadratic programming (QP) problem. Accordingly, the dynamicoptimization included in the OBC may be computationally complex and runon computer servers with general solvers, which may take seconds or evenminutes to produce a solution. Thus, to include OBC on an embeddedsystem for real-time process control, it may be beneficial to improvethe efficiency of OBC while ensuring that it is stabilizing.

Accordingly, one embodiment provides a control method includingdetermining a linear approximation of a pre-determined non-linear modelof a process to be controlled, determining a convex approximation of thenonlinear constraint set, determining an initial stabilizing feasiblecontrol trajectory for a plurality of sample periods of a controltrajectory, executing an optimization-based control algorithm to improvethe initial stabilizing feasible control trajectory for a plurality ofsample periods of a control trajectory, and controlling the controlledprocess by application of the feasible control trajectory within apredetermined time window. In other words, deterministic OBC may beutilized for real-time control of systems with fast dynamics byincluding a stabilization function to produce a stable feasible solution(i.e., a solution that does not increase the cost function) availablefor each predetermined sampling time.

By way of introduction, FIG. 1 depicts an embodiment of a control system10 for a plant/process 12. Generally, the control system 10 may controlthe functioning of the plant/process 12, which may be an industrialmanufacturing system, an automation system, a power generation system, aturbine system, or the like. Accordingly, as depicted, the controlsystem 10 may control the plant/process 12 to transform material inputs14 into material outputs 16. For example, the plant/process 12 may be asugar crystallization process that transforms sugar syrup (i.e.,material input 14) into sugar crystals (i.e., material output 16). Inaddition, the control system 10 may control the output variables (i.e.,controlled variables) 18 by manipulating the input variables 20 (i.e.,manipulated and disturbance variables). Going back to the sugarcrystallization example, the control system 10 may manipulate a steamvalve (i.e., manipulated variable) to control a temperature (i.e.,controlled variable). In some embodiments, the material input can be amanipulated variable as well (for example a controller can control thefeed rate for a material input to the plant).

To optimize the control of the plant/process 12, the control system 10may further include optimization based control (OBC) 22 configured tofind a stabilizing feasible solution for an optimization problem withina predetermined time window. In other words, the OBC 22 may determinefeasible actions (i.e., solution) for the control system 10 to take.Specifically, the OBC 22 may be configured to determine a controltrajectory 26 (i.e., a set of actions) over a control horizon (i.e.,period of time to take the actions). Accordingly, the OBC 22 may samplethe state of the plant/process 12 at specified sampling times. In someembodiments, the state of the plant/process 12 may include the previousoutput variables 18, a desired output trajectory 23, a desired controltrajectory 24, or any combination thereof. Based on the sampled state ofthe plant/process 12, the OBC 22 may determine the control trajectory 26(i.e., a feasible solution to the optimization problem) during thecontrol time. As used herein, control time refers to the time duringwhich the plant/process 12 is functioning, which may be in real-time.After the control trajectory 26 is determined by the OBC 22, in someembodiments, the control trajectory 26 is compared to the desiredcontrol trajectory 24 in a comparator 32 to determine the inputvariables 20 to the plant/process 12 (i.e., actions to be taken in thecontrol system 10). Alternatively, the control trajectory 26 may bedirectly reflected in the input variables 20. It should be appreciatedthat the OBC 22 may be implemented on an embedded system, such asControlLogix, available from available from Rockwell Automation, ofMilwaukee, Wis.

To facilitate determining the control trajectory 26, as depicted, theOBC 22 includes a pre-determined model 28 and a deterministic solver 30.Specifically, the deterministic solver 30 may use a feasible searchstrategy, such as a primal active set method, to determine solutions tothe constrained optimization problem. As will be described in moredetail below, a feasible search strategy begins at a starting pointwithin the feasible region of the control system 10 and moves around thefeasible region to search for an optimum feasible solution (i.e.,control trajectory with minimum cost). In other words, the deterministicsolver 30 may determine various feasible actions (i.e., controltrajectories) that may be taken by the control system 10. Based on thefeasible solutions determined by the deterministic solver 30, the model28 may be utilized to predict the behavior of the process/plant 12. Inlinear systems or non-linear systems with a linear approximation, themodel 28 may be a linear model such as a state space model, a step orimpulse response model, an autoregressive with exogenous terms (ARX)model, a transfer function model, or the like. As such, the OBC 22 maycompare the cost of each feasible solution and select the controltrajectory 26 with the lowest cost.

Ideally, the control trajectory 26 determined is the optimum solutionwith the lowest cost associated, but, as described above, theoptimization calculation may be complex. Accordingly, as will bedescribed in further detail below in the Detailed Example section, thetechniques described herein aim to increase the efficiency of thedynamic optimization calculation. For example, the techniques describedherein may modify an objective (i.e., cost) function to define thecontrol system 10 constraints with simple bounds. Thus, the dynamicoptimization computation may be greatly reduced and executed on anembedded system because many dynamic optimization solvers (e.g.,quadratic-programming (QP) solvers) more efficiently handle simplebounds compared to complex constraints.

Although the dynamic optimization may be efficiently configured, the OBC22 may not always find the optimum (i.e., lowest cost) controltrajectory 26 during each control time. However, in practice, a stablesub-optimal control trajectory 26 may be sufficient. As used herein, thecontrol trajectory 26 is stabilizing when the cost does not increasecompared to the previous step by taking the actions.

To facilitate the functions described herein, it should be appreciatedthat the OBC 22 may include a processor, useful for executing computinginstructions (i.e., steps), and memory, useful for storing computerinstructions (i.e., code) and/or data. As depicted in FIG. 2, the OBC 22may implement the processor through processing circuitry 34 and thememory through memory circuitry 36. More specifically, the processingcircuitry 34 may be configured to handle the general functionality ofthe control system, such as controlling actuators, as well as thefunctions of OBC 22, such as dynamic optimization. In addition, asdepicted in FIG. 3, the processing circuitry 34 may include multipleprocessing components (e.g., parallel processor cores or separateprocessor modules), which may enable the processing circuitry 34 tobetter manage various functions. For example, as depicted, a firstprocessing component 38 may perform the general operations of thecontrol system 10. The general operations of the control system 10 mayinclude controlling components of the control system 10, performingcalculations, and the like. As for the OBC 22 functions, thecomputationally intensive dynamic optimization may be performed on thesecond processing component 40. Accordingly, this enables the dynamicoptimization to be called from the first processing component 38 andexecuted synchronously or asynchronously on the second processingcomponent 40, which may improve the efficiency of the optimizationcalculation. Alternatively, it should be appreciated that the dynamicoptimization may be performed on the first processing core 38 along withthe general functions of the control system 10. Furthermore, asdepicted, the processing circuitry 34 includes N processing components,which may each be configured to handle different functions, such ascalculating a linear approximation, of the control system 10.

Turning back to FIG. 2, the memory circuitry 36 may store computerinstructions (i.e., code) describing the model 28, the deterministicsolver 30, configuration parameters 42, as well as other instructions44, such as computing virtual measurements for unmeasured processvariables for the general functioning of the control system 10.Specifically, the instructions stored in the memory circuit may beconfigured to guide the functioning of the model 28 and thedeterministic solver 30. Accordingly, the memory circuitry 36 iscommunicatively coupled to the processing circuitry 34 to enable toprocessing circuitry 34 to read and/or execute the instructions (i.e.,steps).

Furthermore, the depicted embodiment of the OBC 22 further includes anoutput interface 46, a user interface 48, a network interface 50, and afeedback interface 52. Specifically, the user interface 48 may beconfigured to enable a user to communicate with the OBC 22. For example,as depicted in FIG. 4, the user interface 48 may include agraphical-user-interface (GUI) 54 configured to display metrics of theOBC 22, such as the control trajectory 26 determined. In addition, theuser interface 48 may include buttons 56, which enable the user to inputcommands to the OBC 22. Similar to the user interface 48, the networkinterface 50 may enable a user to communicate with the OBC 22 over anetwork 58, such as a wide-area-network (WAN). In some embodiments, thenetwork 58 may be a EtherNet/IP Network or a ControlNet Network,available from Rockwell Automation, of Milwaukee, Wisc. Morespecifically, as depicted in FIG. 4, the network interface 50 may becommunicatively coupled to the network 58 via a communication module 60.Alternatively, the network interface 50 may be communicatively coupleddirectly the network 58 through the backplane of the OBC 22.Furthermore, as depicted, the network 58 may be communicatively coupledto a remote monitoring/control system 62, such as a supervisory controland data acquisition (SCADA), to enable the user to remotely communicatewith the OBC 22. Accordingly, as depicted in FIG. 2, both the userinterface 48 and the network interface 50 are communicatively coupled tothe processing circuitry 34 to enable user commands to be communicatedto the processing circuitry 34 and information concerning the OBC 22 tobe communicated to the user. Note that each module in memory circuitry36 may be configured such that it can respond as a server responding tothe queries from various interfaces. For example, the model module 28can be queried by the user interface to report its fidelity. In additionthe model module 28 may be called by solver code module 30 to determinethe optimal control trajectory.

Turning back to FIG. 2, as described above, the OBC 22 may be configuredto determine stabilizing feasible control trajectories for the controlsystem 10 based on feedback from the plant/process 12. As such, thefeedback interface 52 may be configured to receive feedback, such as theprevious output variables 18, the desired output trajectory 23, thedesired control trajectory 24, or any combination thereof, andcommunicate it to the processing circuitry 34. For example, the feedbackinterface 52 may be a serial port located on the backplane of the OBC22, which enables the OBC 22 to receive samples from sensors in thecontrol system 10. After the processing circuitry 34 determines acontrol trajectory 26, the control trajectory 26 is communicated to theoutput interface 46. As will be described in more detail below, theprocessing circuitry 34 may utilize various search functions (e.g., QPsolvers) and stabilization functions to determine the control trajectory26. Thus, the output interface 46 may be configured to transmit thecontrol trajectory 26 to the plant/process 12. Similar to the feedbackinterface 52, the output interface 46 may be a serial port located onthe backplane of the OBC 22 to enable the output interface tocommunicate with a controller controlling inputs into the plant/process12. It should be appreciated that as described above, the controller maybe the same processing component, a different core of a processor, or adifferent processor module.

As described above, the processing circuitry 34 may utilize varioussolver methods (i.e., algorithms) to facilitate determining a controltrajectory 26 (i.e., dynamic optimization). Examples of such solvermethods are depicted in FIGS. 5A- 5C. Included in each figure, asdepicted, is a feasible region 64 and an infeasible region 66. Morespecifically, the feasible region 64 is all of the solutions or controltrajectories 26 that do not violate the constraints of the controlsystem 10. On the other hand, the solutions or control trajectories inthe infeasible region 66 violate the constraints of the control system10 and are infeasible. As depicted, constraints are depicted asconstraint lines 68, which separate the feasible region 64 and theinfeasible region 66. As described above, the constraints may be theresult of actuator limits, technological restrictions, economicalrestrictions, and/or safety restrictions.

Specifically, FIGS. 5A and 5B depict feasible search methods (i.e.,algorithms) and FIG. 5C depicts an infeasible search method (i.e.,algorithms). However, it should be appreciated that the figures are notmeant to depict any particular search method or algorithm and are merelyillustrative. As depicted in both FIG. 5A and FIG. 5B, the feasiblesearch method begins from a feasible point 70 within the feasible region64. From the starting feasible point 70, the feasible search methodmoves around the feasible region 64 searching for the optimum solution(i.e., control trajectory) 72 during the control time. In some cases, asshown in FIG. 5A, an optimum control trajectory 72 is found. In others,as shown in FIG. 5B, a suboptimal but still feasible control trajectory74 is found. An example of a feasible search method is a primal activeset solver method (i.e., algorithm). Comparatively, as depicted in FIG.5C, the infeasible search method begins from an infeasible point 76within the infeasible region 66. From the starting infeasible point 76,the infeasible search method determines infeasible solutions until itconverges on the optimum solution (i.e., control trajectory) 72. Anexample of an infeasible search method is a dual active set solvermethod. Accordingly, it should be appreciated that, if terminated beforethe optimum control trajectory 72 is found, the feasible search methodwill produce a feasible control trajectory, but the infeasible searchmethod may produce an infeasible control trajectory.

In addition, as described above, the dynamic optimization (e.g.,feasible search method or infeasible search method) may be runasynchronously from the rest of the control system 10. Thus, if a lessthan optimum control trajectory is found during the control time, theoptimization may continue into the following control times, which givesthe OBC 22 more time for complex optimization calculations. Furthermore,when the optimum control trajectory 72 is found, it may be time shifted,padded, and included in future optimization calculations.

Furthermore, is some cases, the infeasible search method may converge onan optimum solution (i.e., control trajectory) faster than the feasiblesearch method; however, as noted above, the infeasible search mayproduce an infeasible solution if terminated prematurely. Accordingly,to utilize the strengths of both the feasible search method and theinfeasible search method, as depicted in FIG. 6, a sequentialdeterministic optimization based control (OBC) process (i.e., algorithm)78 may be implemented on the OBC 22. Specifically, this may includedividing the maximum execution time (i.e., control time) between theinfeasible search method and the feasible search method and running themsequentially. Beginning the process 78, user-defined setpoints 80 and/orfeasible targets (i.e., infeasible setpoints mapped into feasibleregion) 82 may be input into an initial (i.e., upstream) processing,which may be configured to determine a starting point and constraints(i.e., active set) for the following search methods (process block 84).For example, this may include using a warm start based on a previoussolution. Based on the starting point and constraints, the infeasiblesearch method may run (process block 86). More specifically, theinfeasible search method 86 may run for a period of time shorter thanthe maximum execution time (i.e., control time). If the infeasiblesearch method 86 finds the optimum solution, the process 78 may cease.If on the other hand, during its allotted time, the infeasible search 78has found an infeasible solution, a projection operation occurs toproject the infeasible solution into the feasible region 64 (processblock 88). In some embodiments, this may involve finding the point inthe feasible region 64 closest to the infeasible solution. Anotherexample projection operation 88 may be seen in the Detailed Examplesection.

Because in some embodiments the infeasible search 86 and the projectionoperation 88 may find a suboptimal solution (i.e., control trajectory74) during the allotted execution time, a stabilization process 90 maybe useful, however optional, to stabilize the control system 10. Inother words, the stabilization process 90 may be configured to reducethe risk that the cost of the control system 10 will increase because ofthe suboptimal control trajectory 74. As depicted in FIG. 7, thestabilization process 90 may begin by computing a feasible solution orcontrol trajectory (process block 92). In relation to FIG. 6, thefeasible solution is the projection of the infeasible solution into thefeasible region 64.

To stabilize the control system 10, the feasible solution determined inprocess block 92 (i.e., first stabilizing control trajectory) may becompared to a previous solution (i.e., second stabilizing controltrajectory) determined in process block 94. More specifically, theprevious solution is advanced, which may include shifting and padding,to the next control time (process block 96). As an example, FIG. 8Adepicts a control trajectory 98 with a control horizon of six (i.e., sixtime steps) determined at time i. As described above, the controltrajectory 98 represents a setting of a manipulated variable, such asthe percent a valve is open. Accordingly, based on the controltrajectory 98, the control system 10 will take a first action 100 forthe control time between time i and time i+1. At time i+1, the firstaction 100 has been performed. Thus, as depicted in FIG. 8B, theremaining control trajectory 98 (i.e., trajectory between time i+1 andtime i+6) is time shifted to form a previous solution 102. In addition,because the control horizon for the time shifted control trajectory isfive (i.e., five time steps), the previous solution is 102 is paddedwith an additional control system action 104. In some embodiments, thismay include repeating the last control system action 106. Furthermore,FIG. 8B depicts a newly calculated solution 108, such as by processblock 88, represented by the dashed line. Accordingly, the previoussolution 102 may be compared to the newly calculated solution 108.

Turning back to FIG. 7, after the previous solution 102 and the newlycalculated solution 108 are adjusted to the same control horizon, onecharacteristic to be compared may be the cost of each. Accordingly,based on the objective (i.e. cost) function, the cost for the newlycalculated solution 108 (process block 110) and the cost of the previoussolution 92 may be calculated (process block 112). Next, the cost of thenew solution 110 is compared with the cost of the previous solution 102(process block 114). Finally, the solution (i.e., control trajectory)with the lowest cost may be selected (process block 116).

Turning back to FIG. 6, the feasible search method may be initializedwith either the projected solution determined in process block 88 or thestabilized solution determined in process block 90 (process block 118).Initializing the feasible search 88 with results from the infeasiblesearch 86 may enhance the efficiency of the feasible search 88. Forexample, this may include initializing the feasible search 88 with abetter initial point, initial active set (i.e., set of constraints),and/or matrix factorizations. Furthermore, as will described in moredetail in the Detailed Example section, the searches (i.e., 88 and 86)may be implemented to use the same matrix factorization, which resultsin a seamless transition between the two. Finally, based on the controltrajectory determined by the feasible search 88, the manipulatedvariables 20 may be input into the plant/process 12 to be controlled(process block 120).

FIGS. 9 and 10 depict different embodiments of OBCs 22 utilizing thesequential deterministic optimization based process control 78.Furthermore, as depicted, both embodiments utilize a dual active setsolver method 122 as the infeasible search method and a primal activeset solver method 124 as the feasible search method. As there may beseveral implementations for each solver method (i.e., 122 and 124),different implementations pairing (i.e., an implementation of dualactive search 122 paired with an implementation of primal active setsolver 124) can be found in the Detailed Example section.

As described in FIG. 6, process 78 begins with initial processing 84. Inthe embodiment depicted in FIG. 9, the initial processing 84 may includeinitializing the dual active set solver (process block 126). Asdepicted, initializing the dual active set solver may optionally utilizea warm start 128. Specifically, the warm start 128 may guess an optimalactive set (i.e., set of constraints) 130 based on a previous solution132. Accordingly, the warm start may be useful when the dynamicoptimization problem is similar in subsequent control times. If a warmstart is not utilized, the initial active set is empty. Based on theinitial active set, a dual feasible point and active set 134 may bedetermined. If a warm start 128 is used, constraints corresponding todual variables (i.e., Lagrange multipliers) with the wrong (i.e.,negative) sign in the guessed optimal active set 130 may be removeduntil all dual variable have correct (i.e., nonnegative) sign. A furtherdescription of the dual variables may be seen in the Detailed Examplesection. Alternatively, if the warm start 128 is not used, the dualfeasible point may be the unconstrained solution of the optimizationproblem.

Next, as described in FIG. 6, the dual feasible point and active set 134are used to initialize the dual active set solver method 122, and if thedual active set solver method 122 finds the optimum solution in itsallotted time, the process 78 terminates. If the allotted time hasexpired, the projection processing 88 projects the infeasible solutioninto the feasible region 64. In the depicted embodiment, the projectionprocessing produces a primal feasible point and active set 136.

The primal active set solver method 124 may then be initialized throughvarious manners. For example, the primal active set solver method 124may use the result of the projection processing 136, adjust the resultof the projection processing 136 by adding selected constraints thatbecome active after the projection operation 136, or adjust the resultof the projection processing 136 by adding selected constraints thatbecome active after the projection operation 136 and by addingartificial constraints to satisfy a constrained sub-problem. Finally,the primal active set solver method 124 may use the remaining time tosearch for an optimum solution (i.e., control trajectory).

In addition to the functions included in FIG. 9, FIG. 10 furtherincludes the stabilization process 90, which may be useful foroptimization problems with simple bound constraints. Specially, for thistype of well-structured constraints, the projection operation can beefficiently calculated, for example by a clipping operation.

Furthermore, as described above, the OBC 22 may be run on multiplecomponents (i.e., 38 and 40), which enables multiple processes to be runasynchronously. Accordingly, as depicted in FIG. 11, the OBC 22 mayutilize multiple optimization searches in parallel, which as should beappreciate may increase the probability of finding the optimum solution(i.e. control trajectory). Specifically, FIG. 11 depicts a paralleldeterministic optimization based control (OBC) process 138 including thesequential deterministic optimization based control (OBC) process 78. Asdepicted, the sequential deterministic optimization based control (OBC)process 78 is in parallel with another optimization search (i.e.,solver) 140, which may be any optimization solver method (i.e.,algorithm). It may be useful to parallel process 78 with anotherinfeasible search 86 that is allotted a longer maximum execution timebecause, as described above, process 78 will return a feasible solution.Alternatively, process 78 may be paralleled with another feasible search88 so that the feasible search may serve as a backup to process 78 whileadditionally searching for the optimum solution (i.e., controltrajectory).

After process 78 and optimization search 140 have terminated, theresults may be compared. For example, as depicted, the cost of process78 (process block 142) and the cost of the parallel optimization search140 (process block 144) may be calculated based on the objective (i.e.,cost) function. These costs may be compared (process block 146) and thelower cost may be selected (process block 148). Finally, based on thecontrol trajectory selected, the manipulated variables 20 input into theplant/process 12 may be controlled (process block 120).

In addition, it should be appreciated that the optimization searches(i.e., 78 and 140) may be run synchronously or asynchronously. Forexample, process 78 and parallel optimization search 140 may be calledsubstantially simultaneously. Alternatively, one (i.e., 78 or 140) maybe called and allowed to run on a processing component (e.g., 38 or 40)while the control system 10 continues with its normal functioning, whichprovides the search method (i.e., 78 or 140) with a longer search time.Then, the other (i.e., 78 or 140) may be called subsequently.

DETAILED EXAMPLE

Below is a detailed example to help illustrate the techniques taughtherein. First, as described above, the efficiency of the OBC 22 may beincreased by simplifying the dynamic optimization calculation. In someembodiments, this may include simplifying the model 28. Generally, themodel 28 may be expressed as a standard formula.

y=A _(u) Δu+y _(free)  (1)

where

Δu=[(Δu₁)^(T), (Δu₂)^(T), . . . , (Δu_(n) _(u) )^(T)]^(T)—future changesof all inputs arranged in a vector

Δu_(n)=[Δu_(n)(t), Δu_(n)(t+1), . . . , Δu_(n)(t+n_(c)−1)]^(T)—futurechanges of n-th input

y=[(y₁)^(T),(y₂)^(T), . . . , (y_(n) _(y) )^(T)]^(T)—all predictedoutput trajectories arranged in a vector

y_(j)=[y_(j)(t+1), . . . , y_(j)(t+p)]^(T)—predicted trajectory of j-thoutput)

y_(free)=[(y¹ _(free))^(T), (y² _(free))^(T), . . . , (y^(n) ^(y)_(free))^(T)]^(T)—all free response trajectories arranged in a vector

y^(j) _(free)=[y^(j) _(free)(t+1), . . . , y^(j) _(free)(t+p)]^(T)—freeresponse trajectory of j-th output

n_(u)—number of inputs

n_(y)—number of outputs

n_(c)—control horizon

p—prediction horizont

N=n_(c)—number of control moves

Specifically, the standard formula (equation 1) reflects thesuperposition principle: the total future response (i.e., outputtrajectory) is a sum of responses to future inputs and responses to pastinputs, disturbances, and initial conditions (i.e., y_(free)). Theresponse y_(free) may be updated at each sample time according toparticular model structure. Accordingly, the predicted output trajectoryof the j-th output may be expressed similar to equation 1.

y _(j) =A _(uj) Δu+y ^(j) _(free)  (2)

where

A_(uj)=[A_(uj) ¹, A_(uj) ², . . . , A_(uj) ^(n) ^(u) ]—prediction matrixof j-th output

A_(u)=[A_(u1) ^(T), A_(u2) ^(T), . . . , A_(un) _(y) ^(T),]^(T)—prediction matrix

Furthermore, the following equation represents the changes of inputs(i.e., incremental inputs).

Δu(k)=u(k)−u(k−1)  (3)

In addition, the following represent the absolute inputs.

u _(n) =[u _(n)(t),u _(n)(t+1), . . . , u _(n)(t+n _(c)−1)]^(T)—futuretrajectory of n-th input

u=[(u ₁)^(T), (u ₂)^(T), . . . , (u _(n) _(u) )]^(T)—all future inputtrajectories arranged in a vector  (4)

Thus, to help simplify the above equations, relative inputs may be usedinstead of incremental inputs. Specifically, the constraints may besimplified and the Hessian matrix may be better conditioned. As usedherein u_(r) represents the relative input and is defined as differencefrom current value of input. Accordingly, u_(r) may be expressed asfollows.

u _(ri)(k)=u _(i)(k)−u _(0i)  (5)

where

u_(ri)—i-th relative input

u_(0i)=u_(i)(t−1)—current value of the i-th input

Based on the definition of u_(n), equation (5) may alternativelyexpressed as follows.

u _(ri) =u _(i)−1u _(0i), 1^(T)=[1,1, . . . 1]^(T)  (6)

Thus, the linear transformation between the relative input (u_(r)) andthe incremental input (Δu) is as follows.

$\begin{matrix}{{{{\Delta {u_{i}(k)}} = {{{u_{i}(k)} - {u_{i}\left( {k - 1} \right)}} = {{\Delta {u_{ri}(k)}} = {{u_{ri}(k)} - {u_{ri}\left( {k - 1} \right)}}}}}{{u_{ri}(k)} = {\sum\limits_{j = t}^{k}{\Delta {u_{i}(j)}}}}}{{u_{ri} = {S\Delta u_{i}}},{S = \begin{bmatrix}1 & 0 & 0 & \ldots & 0 \\1 & 1 & 0 & \; & 0 \\1 & 1 & 1 & \; & 0 \\\vdots & \; & \; & \ddots & 0 \\1 & 1 & 1 & 1 & 1\end{bmatrix}},{{\Delta \; u_{i}} = {S^{- 1}u_{ri}}},{S^{- 1} = \begin{bmatrix}1 & 0 & 0 & \ldots & 0 \\{- 1} & 1 & 0 & \; & 0 \\0 & {- 1} & 1 & \; & \; \\\vdots & \; & \; & \ddots & \vdots \\0 & 0 & \ldots & {- 1} & 1\end{bmatrix}}}{{{\Delta \; u} = {S_{c}^{- 1}u_{r}}},{S_{c}^{- 1} = \begin{bmatrix}S^{- 1} & 0 & \ldots & 0 \\0 & S^{- 1} & \; & 0 \\\vdots & \; & \ddots & \; \\0 & 0 & \; & S^{- 1}\end{bmatrix}},{u_{r} = {S_{c}\Delta u}},{S_{c} = \begin{bmatrix}S & 0 & \ldots & 0 \\0 & S & \; & 0 \\\vdots & \; & \ddots & \; \\0 & 0 & \; & S\end{bmatrix}}}{{u_{r} = {u - {I_{c}u_{0}}}},{I_{c} = \begin{bmatrix}1 & 0 & \ldots & 0 \\0 & 1 & \; & 0 \\\vdots & \; & \ddots & \; \\0 & 0 & \; & 1\end{bmatrix}},{u_{0} = \begin{bmatrix}u_{01} \\u_{02} \\\vdots \\u_{0n_{u}}\end{bmatrix}}}{u = {{S_{c}{\Delta u}} + {I_{c}u_{0}}}}} & (7)\end{matrix}$

Accordingly, the model 28 (i.e., equation 1) based on the relative inputmay be expressed as follows.

y=A _(u) Δu+y _(free)=(A _(u) S _(c) ⁻¹)u _(r) +y _(free) =A _(r) u _(r)+y _(free)  (8)

In addition to simplifying the model 28, cost function (i.e., objectionfunction) may also be simplified to enhance the efficiency of the OBC22. Generally, a quadratic objective function may be expressed asfollows.

$\begin{matrix}{J_{y} = {{\sum\limits_{j = 1}^{n_{y}}{\left( {y_{j} - y_{j}^{t}} \right)^{T}{W_{j}\left( {y_{j} - y_{j}^{t}} \right)}}} + {\sum\limits_{n = 1}^{n_{u}}{\Delta \; u_{n}^{T}B_{n}\Delta \; u_{n}}} + {\sum\limits_{n = 1}^{n_{u}}{\left( {u_{n} - u_{n}^{t}} \right)^{T}{V_{n}\left( {u_{n} - u_{n}^{t}} \right)}}}}} & (9)\end{matrix}$

where

W_(j), B_(n), V_(n)-positive semidefinite weight matrices (tuningparameters)

y^(t) _(j)=[y^(t) _(j)(t+1), y^(t) _(j)(t+2), . . . , y^(t)_(j)(t+p)]^(T)—output target trajectory

u_(n) ^(t)=[u_(n) ^(t)(t+1), u_(n) ^(t)(t+2), . . . , u_(n)^(t)(t+p)]^(T)—input target trajectory

Using standard algebraic manipulation the objective function can betransformed to the following quadratic form.

J _(y)(Δu)=½Δu ^(T) G _(y) Δu+f _(y) ^(T) Δu+C  (10)

where G_(y) is Hessian, f_(y) is a linear term and C is constant termthat can be omitted for optimization. We assume that B_(n) matrices arepositive definite and so G_(y) is also positive definite.

One of the most important features of OBC is that it can handleconstraints. The constraints for the n-th input may be represented asfollows.

i=0 . . . n _(c)−1:

u _(min)(n)≤u _(n)(t+i)≤u _(max)(n)—simple bounds

Δu _(min)(n)≤Δu _(n)(t+i)≤Δu _(max)(n)—rate of change constraints  (11)

In a compact form, the constraints may be represented as follows.

Δu _(max) =[Δu _(max)(1), . . . , Δu _(max)(n _(u))]^(T)

Δu _(min) =[Δu _(min)(1), . . . , Δu _(min)(n _(u))]^(T)

I_(c)Δu_(min)≤Δu≤I_(c)Δu_(max)

u _(max) =[u _(max)(1), . . . , u _(max)(n _(u))]^(T)

u _(min) =[u _(min)(1), . . . , u _(min)(n _(u))]^(T)

u_(min)≤u≤u_(max)

u _(min) ≤S _(c) Δu+I _(c) u ₀ ≤u _(max)

I _(c)(u _(min) −u ₀)≤S _(c) Δu≤I _(c)(u _(max) −u ₀)  (12)

In some cases, controlled variable (i.e., output) may lie within aspecific range (i.e., between y_(min) and y_(max)). This formulation maylead to more robust and less aggressive control actions but may alsoresults in a more complex optimization calculation. However, if theoptimization calculation has hard constraints on outputs, then aninfeasible control trajectory may be found. Accordingly, softconstraints with slack variables may be used. However, in order to notto dramatically increase number of decision variables or number ofcomplex constraints, the optimization calculation may be simplified byonly adding one additional variable with simple bounds per constrainedoutput. This additional variable may be interpreted as an optimizedtarget that must lie within the specific output range (i.e., betweeny_(min) and y_(max)). For example, output target trajectory) y^(t) maybe defined by an additional decision variable z and the constrainedoutputs may be defined by n_(z) with indices k₁ . . . k_(nz).Accordingly, objective function and constraints may be defined asfollows.

$\begin{matrix}{{{J_{z}\left( {{\Delta \; u},z} \right)} = {{\sum\limits_{i = 1}^{n_{z}}{\left( {y_{k_{i}} - {1z_{i}}} \right)^{T}{M_{i}\left( {y_{k_{i}} - {1z_{i}}} \right)}}} = {{\begin{bmatrix}{\Delta \; u} \\z\end{bmatrix}^{T}{G_{z}\begin{bmatrix}{\Delta \; u} \\z\end{bmatrix}}} + {f_{z}^{T}\begin{bmatrix}{\Delta \; u} \\z\end{bmatrix}} + C}}}\mspace{20mu} {y_{\min} \leq z \leq y_{\max}}\mspace{20mu} {{1 = \left\lbrack {1,1,{\ldots \; 1}} \right\rbrack^{T}},{z = \left\lbrack {z_{1},z_{2},{\ldots \; z_{n_{z}}}} \right\rbrack^{T}}}{{y_{\max} = \left\lbrack {{y_{\max}(1)},\ldots \;,{y_{\max}\left( n_{z} \right)}} \right\rbrack^{T}},{y_{\min} = \left\lbrack {{y_{\min}(1)},\ldots \;,{y_{\min}\left( n_{z} \right)}} \right\rbrack^{T}}}} & (13)\end{matrix}$

Furthermore, because overall objective function J is then sumJ=J_(y)+J_(z) the optimization calculation may be as follows.

$\begin{matrix}{{{\min\limits_{{\Delta \; u},z}\; {J\left( {{\Delta \; u},z} \right)}} = {{\begin{bmatrix}{\Delta \; u} \\z\end{bmatrix}^{T}{G\begin{bmatrix}{\Delta \; u} \\z\end{bmatrix}}} + {f^{T}\begin{bmatrix}{\Delta \; u} \\z\end{bmatrix}} + C}}{{I_{c}\Delta \; u_{\min}} \leq {\Delta \; u} \leq {I_{c}\Delta \; u_{\max}}}{{I_{c}\left( {u_{\min} - u_{0}} \right)} \leq {S_{c}\Delta \; u} \leq {I_{c}\left( {u_{\max} - u_{0}} \right)}}{z_{\min} \leq z \leq z_{\max}}} & (14)\end{matrix}$

Using the transformation Δu=S_(c) ⁻¹u_(r) in equation (7), theoptimization problem based on relative inputs (u_(r)) is as follows.

$\begin{matrix}{{{\min\limits_{u_{r},z}\; {J_{r}\left( {u_{r},z} \right)}} = {{\begin{bmatrix}u_{r} \\z\end{bmatrix}^{T}{G_{r}\begin{bmatrix}u_{r} \\z\end{bmatrix}}} + {f_{r}^{T}\begin{bmatrix}u_{r} \\z\end{bmatrix}}}}{{I_{c}\left( {u_{\min} - u_{0}} \right)} \leq u_{r} \leq {I_{c}\left( \; {u_{\max} - u_{0}} \right)}}{{I_{c}\Delta \; u_{\min}} \leq {S_{c}^{- 1}\; u_{r}} \leq {I_{c}\Delta \; u_{\max}}}{z_{\min} \leq z \leq z_{\max}}} & (15)\end{matrix}$

Accordingly, taking advantage of that fact that feasible search methodmore efficiently hand simple bound constraints as compared to complexconstraints 150 (e.g., a rate of change constraint), depicted in FIG.12A, the simple bound constraints 152 may be configured to act as afunnel as depicted in FIG. 12B. Thus, the simple bounds may be expressedas follows.

$\begin{matrix}{\mspace{79mu} {{{n = {{1\ldots \; {n_{u}:\mspace{76mu} {u_{\max}^{rn}(i)}}} = {\min \left( {{{u_{\max}(n)} - {u_{0}(n)}},{{i \cdot \Delta}\; {u_{\max}(n)}}} \right)}}},{i = {1\mspace{11mu} \ldots \mspace{11mu} n_{c}}}}\mspace{76mu} \underset{\_}{{u_{\min}^{rn}(i)} = {\max \left( {{{u_{\min}(n)} - {u_{0}(n)}},{{{- i} \cdot \Delta}\; {u_{\min}(n)}}} \right)}}{{u_{\max}^{rn} = \left\lbrack {{u_{\max}^{rn}(1)},\ldots \;,{u_{\max}^{rn}\left( n_{C} \right)}} \right\rbrack^{T}},{u_{\max}^{r} = \left\lbrack {\left( u_{\max}^{r1} \right)^{T},\ldots \;,\left( u_{\max}^{{rn}_{u}} \right)^{T}} \right\rbrack^{T}}}{{u_{\min}^{rn} = \left\lbrack {{u_{\max}^{rn}(1)},\ldots \;,{u_{\max}^{rn}\left( n_{C} \right)}} \right\rbrack^{T}},{u_{\min}^{r} = \left\lbrack {\left( u_{\min}^{r1} \right)^{T},\ldots \;,\left( u_{\min}^{{rn}_{u}} \right)^{T}} \right\rbrack^{T}}}}} & (16)\end{matrix}$

Thus, the optimization calculation may be expressed as follows.

$\begin{matrix}{{{\min\limits_{u_{r},z}{J_{r}\left( {u_{r},z} \right)}} = {{\begin{bmatrix}u_{r} \\z\end{bmatrix}^{T}{G_{r}\begin{bmatrix}u_{r} \\z\end{bmatrix}}} + {f_{r}^{T}\begin{bmatrix}u_{r} \\z\end{bmatrix}}}}{u_{\min}^{r} \leq u_{r} \leq u_{\max}^{r}}{{{I_{c}{\Delta u}_{\min}} \leq {Ru}_{r} \leq {I_{c}{\Delta u}_{\max}}},\text{}{z_{\min} \leq z \leq z_{\max}}}{{R = \begin{bmatrix}R_{u} & 0 & \cdots & 0 \\0 & R_{u} & \; & 0 \\\vdots & \; & \ddots & \; \\0 & 0 & \; & R_{u}\end{bmatrix}},{R_{u} = \begin{bmatrix}{- 1} & 1 & 0 & \; & 0 \\0 & {- 1} & 1 & \; & 0 \\\; & \; & \; & \ddots & \vdots \\0 & 0 & \cdots & {- 1} & 1\end{bmatrix}},}} & (17)\end{matrix}$

Furthermore, using a common variable x=[u_(r) ^(T), z^(T)]^(T) with nelements we get the following resulting optimization problem in astandard form as follows.

$\begin{matrix}{{{{\min\limits_{x}{J_{r}(x)}} = {{x^{T}G_{r}x} + {f_{r}^{T}x}}}{x_{{mi}n} \leq x \leq x_{\max}}{{\Delta x}_{\min} \leq {Rx} \leq {\Delta x}_{\max}}{R = \begin{bmatrix}R_{u} & 0 & \cdots & 0 \\0 & R_{u} & \; & 0 \\\vdots & \; & \ddots & \; \\0 & 0 & \; & R_{u}\end{bmatrix}},{R:{{n_{u}\left( {N - 1} \right)} \times n}},{R_{u} = \begin{bmatrix}{- 1} & 1 & 0 & \; & 0 & \cdots & 0 \\0 & {- 1} & 1 & \; & 0 & \; & 0 \\\; & \; & \; & \ddots & \vdots & \; & \; \\0 & 0 & \cdots & {- 1} & 1 & \; & 0\end{bmatrix}},{R_{u}:{\left( {N - 1} \right) \times n}}}{n = {{n_{u}N} + n_{z}}}} & (18)\end{matrix}$

where

n—number of decision variables (length of x) p n_(u)—number of inputs(manipulated variables)

N—number of control moves (blocks) per control input

n_(z)—number of constrained outputs (controlled variables)

R—delta constraints

Further efficiency enhancements may be made to the various searchmethod. For example, rotations may be used to simplify an active set(i.e., constraints) matrix A. Specifically, factorization may be used torotate the active set matrix when a constraint is added or deleted fromactive set. The described stable implementation is based on a modifiedCholesky factorization of the Hessian matrix G where U is an uppertriangular matrix.

G=UU^(T)  (19)

Thus, the modified Cholesky factorization may be computed using astandard Cholesky factorization and a permutation matrix P as follows.

$\begin{matrix}{{{{PGP^{T}} = {L_{2}L_{2}^{T}}}{G = {UU^{T}}}},{U = {P^{T}L_{2}P}},{P = \begin{bmatrix}\; & \; & \; & ⋰ \\0 & 0 & 1 & \; \\0 & 1 & 0 & \; \\1 & 0 & 0 & \;\end{bmatrix}}} & (20)\end{matrix}$

As described in equation (20) the Hessian matrix is first permuted. Thenthe standard Cholesky factorization is calculated. Finally, thetriangular matrix is permuted to get matrix U. Accordingly the followingfactorization of the matrix U⁻¹N⁺, where T⁺ is an upper triangularmatrix, may be used.

$\begin{matrix}{{{\left( M^{+} \right)^{T}U^{1}N^{+}} = \begin{bmatrix}T^{+} \\0\end{bmatrix}}{{\left( M^{+} \right)^{T}M^{+}} = D^{- 1}}{D = {{{diag}\left( {d_{1},d_{2},{..d_{n}}} \right)} = {\begin{bmatrix}D_{l} & 0 \\0 & D_{2}\end{bmatrix}\begin{matrix}{\} q} \\{{\} n} - q}\end{matrix}}}}} & (21)\end{matrix}$

Furthermore, the J⁺ matrix may be calculated by the following equation.

$\begin{matrix}{J^{+} = {{U^{- T}M^{+}} = {\left\lbrack {L^{- T}M_{1}^{+}\ L^{- T}M_{2}^{+}} \right\rbrack = \left\lbrack {\underset{q}{\underset{}{J_{1}^{+}}}\underset{n - q}{\underset{}{J_{2}^{+}}}} \right\rbrack}}} & (22)\end{matrix}$

The described factorization has advantage of upper triangularity ofmatrix U⁻¹. Thus when N⁺ is composed of an initial part of an identitymatrix the factorization is readily available.

As described above, plane rotations may be performed when a constraintis added or removed from the active set. For example, constraint p maybe removed from the active set because the corresponding dual variableu_(p) is negative. Furthermore, the following matrix structure isassumed.

N ⁺ =[N ₁ n _(p) N ₂ ], N=[N ₁ N ₂],   (23)

Matrix N⁺ may first be reordered so that normal n_(p) is moved to thelast column as follows.

N _(r) ⁺ =[N ₁ N ₂ n _(p)]  (24)

Thus, the triangularity of the matrix T⁺ may be affected as follows.

$\begin{matrix}{{\left( J^{+} \right)^{T}N_{r}^{+}} = {{\begin{bmatrix}T_{H}^{+} \\0\end{bmatrix}\begin{matrix}{\} q} \\{{\} n} - q}\end{matrix}} = {\begin{bmatrix}T_{1} & S & t_{p1} \\0 & V & t_{p2} \\0 & 0 & 0\end{bmatrix}\begin{matrix}{{\} p} - 1} \\{{\} q} - p + 1} \\{{\} n} - q}\end{matrix}}}} & (25)\end{matrix}$

In order to restore its upper triangularity, a set of plane rotations Qmay be applied to upper Hessenberg sub-matrix [V t_(p2)] where γ isscalar.

$\begin{matrix}{{{Q\left\lbrack {V\ t_{p2}} \right\rbrack} = \begin{bmatrix}T_{2} & d_{1} \\0 & J^{1}\end{bmatrix}},{T_{r}^{+} = {{\begin{bmatrix}1 & 0 \\0 & Q\end{bmatrix}T_{H}^{+}} = {\begin{bmatrix}T_{1} & S & d_{1} \\0 & T_{2} & \; \\0 & 0 & \gamma\end{bmatrix} = \begin{bmatrix}T & d_{1} \\0 & \gamma\end{bmatrix}}}},{d = \begin{bmatrix}d_{1} \\\gamma\end{bmatrix}}} & (26)\end{matrix}$

To nullify elements under the main diagonal, the rotations may beapplied to matrix rows in the following order: (1,2), (2,3), . . . (q−p,q−p+1). The same set of plane rotations Q may be applied tocorresponding columns (p,p+1), . . . (q−1,q) of matrix J⁺ as follows.

$\begin{matrix}{{J^{+} = {\left\lbrack {\underset{q}{\underset{}{J_{1}^{+}}}\underset{n - q}{\underset{}{J_{2}^{+}}}} \right\rbrack = \left\lbrack {\underset{p}{\underset{}{J_{1}^{+}}}\underset{q - p}{\underset{}{J_{12}^{+}}}\underset{n - q}{\underset{}{J_{2}^{+}}}} \right\rbrack}}{J = {\left\lbrack {\underset{q}{\underset{}{J_{11}^{+}J_{12}^{+}Q^{T}}}\underset{n - q}{\underset{}{J_{2}^{+}}}} \right\rbrack = \left\lbrack {\underset{q - 1}{\underset{}{J_{1}}}\underset{n - q + 1}{\underset{}{J_{2}}}} \right\rbrack}}} & (27)\end{matrix}$

On the other hand, constraints may be added to the active set when. Forexample, constraint Ilk may be added to active set as follows.

A ₂ =AU{k}, N ₂ =[N n _(k) ], A ⁺ =A ₂ U{p}, N ⁺ =[N ₂ n _(p)]  (28)

To update the matrix T_(r) ⁺, the last column d is shifted to the rightand a new column h is put in its place as follows.

$\begin{matrix}{{N^{+} = \left\lbrack {N\mspace{14mu} n_{k}\mspace{14mu} n_{p}} \right\rbrack}{u_{p} = {\left( {1 - t} \right)u_{p}}}{T_{r}^{+} = \left\lbrack {{T_{r}^{+}\left( {:{,{1:q}}} \right)},\ h,\ d} \right\rbrack}{{h = {{J^{T}n_{k}} = {\begin{bmatrix}h_{1} \\h_{2}\end{bmatrix}\begin{matrix}{\} q} \\{{\} n} - q}\end{matrix}}}},{d = {\begin{bmatrix}d_{1} \\d_{2}\end{bmatrix}\begin{matrix}{\} q} \\{{\} n} - q}\end{matrix}}},{d_{2} = \left( \underset{n - q}{\underset{}{\left\lbrack {\gamma \mspace{14mu} 0{\ldots 0}} \right\rbrack}} \right)^{T}}}} & (29)\end{matrix}$

Similar to removing a constraint, for the new active set, a sequence ofrotations Q may be applied to restore the triangularity. Based on thestructure of the T_(r) ⁺, it may be sufficient to nullify all but thefirst element of the vector h₂. Specifically, the rotations to beapplied in the planes (n−q, n−q−1), (n−q−1, n−q−2), . . . (2, 1) asfollows.

$\begin{matrix}{{h_{2} = {\left. \begin{bmatrix}x \\x \\x \\x\end{bmatrix}\rightarrow\left. \begin{bmatrix}x \\x \\y \\0\end{bmatrix}\rightarrow\left. \begin{bmatrix}x \\z \\0 \\0\end{bmatrix}\rightarrow\begin{bmatrix}\alpha \\0 \\0 \\0\end{bmatrix} \right. \right. \right. = \overset{¯}{h_{2}}}}{d_{2} = {\left. \begin{bmatrix}x \\0 \\0 \\0\end{bmatrix}\rightarrow\left. \begin{bmatrix}x \\0 \\0 \\0\end{bmatrix}\rightarrow\left. \begin{bmatrix}x \\0 \\0 \\0\end{bmatrix}\rightarrow\begin{bmatrix}a \\b \\0 \\0\end{bmatrix} \right. \right. \right. = \overset{¯}{d_{2}}}}} & (30)\end{matrix}$

After these rotations are applied to T_(r) ⁺, the resulting triangularmatrix is denoted as T⁺. If fast rotations are applied then alsocorresponding diagonal elements of matrix D are updated. The same set ofrotations may be applied to the columns of matrix J⁺ as follows.

$\begin{matrix}{{{{{\left. {{Qh_{2}} = \begin{bmatrix}\alpha \\0\end{bmatrix}} \right\} n} - q} = \overset{¯}{h_{2}}},{{Qd_{2}} = \overset{¯}{d_{2}}}}{J = \left\lbrack {\underset{q}{\underset{}{J_{1}^{+}}}\underset{n - q}{\underset{}{J_{2}^{+}Q^{T}}}} \right\rbrack}} & (31)\end{matrix}$

After the above operations, the number of active constraints andnotation of matrices may be updated as follows.

$\begin{matrix}{{q = {q + 1}},{J^{+} = {J = \left\lbrack {\underset{q}{\underset{}{J_{1}^{+}}}\; \underset{n - q}{\underset{}{J_{2}^{+}}}} \right\rbrack}},} & (32)\end{matrix}$

In a dual active set solver method, when constraint n⁺ is violated, thestep direction may be calculated as follows.

$\begin{matrix}{{{J^{T}n^{+}} = {{\begin{bmatrix}J_{1}^{T} \\J_{2}^{T}\end{bmatrix}n^{+}} = {\begin{bmatrix}d_{1} \\d_{2}\end{bmatrix} = d}}}{z = {{Hn^{+}} = {J_{2}D_{2}d_{2}}}}{r = {{- T^{1}}d_{1}}}} & (33)\end{matrix}$

As used herein the primal space is represented by z and the dual spaceis represented by r. And, when the constraint n⁺ is selected to bedeleted from the active set, step direction may be calculated asfollows.

z=−Hn ⁺ =−J ₂ D ₂ d ₂

r=T⁻¹d₁  (34)

As described in the Detailed Description section, the dual active setsolver 122 and the primal active set solver 124 search methods may beimplemented through various implementations. More specifically, eachimplementation of the solver methods (i.e., 122 and 124) may be based ona method of matrix factorization. For example, the techniques describedherein may be based on the Goldfarb method of matrix factorization.Alternatively, other implementations of the search methods may be used;however, to improve the transition between the dual active set solver122 and the primal active set solver 124, it may be beneficial toutilize implementations that are based on the same matrix factorizationmethod. Thus, the matrix factorization function may be shared by thesolver methods (i.e., 122 and 124) and the transition may be seamless.Other matrix factorization pairings that facilitate these benefits areas follows:

Dual: Goldfarb and Idnani method, Primal: Goldfarb method;

Dual: QPSchur method, Primal: Schur complement method;

Dual: Gram-Schmidt range-space method, Primal: Gram-Schmidt range-spacemethod; or

Dual: Null-space method, Primal: Null-space method.

Finally, as described above in the Detailed Description Section, variousembodiments of the projection operation 136 may be used. For example,when constraints are simple, sparse or well-structured the projectionoperation 136 may be efficiently calculated. In particular, theprojection operation (PROJ) 136 described below may be used on thefollowing simple constraints.

i=0 . . . n _(c)−1:

u _(min)(n)≤u _(n)(t+i)≤u _(max)(n)—simple bounds

Δu _(min)(n)≤Δu _(n)(t+i)≤Δu _(max)(n)—rate of change constraints  (35)

Hard constraints may be imposed by clipping as follows.

$\begin{matrix}{{C_{abs}\left( {u_{n}\left( {t + i} \right)} \right)} = \left\{ \begin{matrix}{u_{\min}(n)} & {{{if}\mspace{14mu} {u_{n}\left( {t + i} \right)}} < {u_{\min}(n)}} \\{u_{\max}(n)} & {{{if}\mspace{14mu} {u_{\max}\left( {t + i} \right)}} < {u_{n}\left( {t + i} \right)}} \\{u_{n}\left( {t + i} \right)} & {otherwise}\end{matrix} \right.} & (36)\end{matrix}$

Similarly, rate of change (ROC) constraints may be imposed by clippingas follows

$\begin{matrix}{{C_{roc}\left( {u_{n}\left( {t + i} \right)} \right)} = \left\{ \begin{matrix}{{u_{n}(t)} + {{\Delta u}_{\min}(n)}} & {{{{if}\mspace{14mu} {u_{n}\left( {t + i} \right)}} - {u_{n}(t)}} < {{\Delta u}_{\min}(n)}} \\{{u_{n}(t)} + {{\Delta u}_{\max}(n)}} & {{{{if}\mspace{14mu} {u_{n}\left( {t + i} \right)}} - {u_{n}(t)}} > {{\Delta u}_{\max}(n)}} \\{u_{n}\left( {t + i} \right)} & {otherwise}\end{matrix} \right.} & (37)\end{matrix}$

Thus, the projection operation 136 may be defined as follows.

PROJ(u _(n))=C _(roc)(C _(abs)(u _(n)))  (38)

Once a feasible point is obtained, the corresponding initial active setmay be calculated. Furthermore, to exploit the final matrixfactorizations of infeasible search 86, the initial active set of thefeasible solver 88 may be equal to the final active set of infeasiblesearch 86. After the projection operation 136, some of the constraintsmay no longer be active and vice versa. If a constraint is no longeractive, then it may be designated as artificial. During the course ofthe run of the optimization calculation, artificial active constraintsmay be deleted preferentially. If a constraint becomes active after theprojection operation, it is not added to the initial active set of thefeasible solver 88, but may be added during the execution of thefeasible solver 88.

Generally, the above techniques enable deterministic optimizationcontrol to be used with plants/process 12 with fast dynamics. Morespecifically, the OBC 22 is configured to advantageous utilize thecharacteristics of feasible and infeasible search methods to provide afeasible control trajectory during each control time, which in someembodiments includes using the infeasible search method and the feasiblesearch method sequentially. In other words, the control time may bedivided between the two search methods.

While only certain features of the invention have been illustrated anddescribed herein, many modifications and changes will occur to thoseskilled in the art. It is, therefore, to be understood that the appendedclaims are intended to cover all such modifications and changes as fallwithin the true spirit of the invention.

1.-20. (canceled)
 21. A method of controlling operation of an industrialautomation process, comprising: performing, using a control system, afirst search during a predetermined sample period based at least in parton one or more sensor samples indicative of a state of the industrialautomation process to determine a first control trajectory of aplurality of input variables of the industrial automation process;performing, using the control system, a second search during thepredetermined sample period based at least in part on the one or moresensor samples indicative of the state of the industrial automationprocess to determine a second control trajectory of the input variablesof the industrial automation process; selecting the first controltrajectory for instructing control actions when the first controltrajectory is determined within a desired time and is an optimal controltrajectory; and selecting the second control trajectory for instructingcontrol actions when the first control trajectory is not selected. 22.The method of claim 21, wherein the first search comprises an infeasiblesearch beginning at a first point outside a feasible region for controltrajectories.
 23. The method of claim 21, wherein the second searchcomprises a feasible search beginning at a second point inside thefeasible region for control trajectories.
 24. The method of claim 21,wherein the first search converges on an optimal control trajectoryfaster than the second search.
 25. The method of claim 21, wherein thefirst and second searches are run sequentially in the sample period. 26.The method of claim 21, wherein the first and second searches areperformed starting from a solution from a previous sample period. 27.The method of claim 21, wherein the first and second searches areperformed asynchronously from other operations of the control system.28. The method of claim 21, wherein: performing the first searchcomprises solving a constrained quadratic programming problem of anoptimization-based control algorithm; and performing the second searchcomprises solving the constrained quadratic programming problem.
 29. Themethod of claim 21, wherein: performing the first search comprisesexecuting a dual active set solver algorithm to solve a constrainedquadratic programming problem; and performing the second searchcomprises executing a primal active set solver algorithm to solve theconstrained quadratic programming problem.
 30. A system of controllingan industrial automation process, comprising: memory configured to storeexecutable code; and processing circuitry communicatively coupled to thememory, wherein the processing circuitry is configured to execute theexecutable code to: perform, using a control system, a first searchduring a predetermined sample period based at least in part on one ormore sensor samples indicative of a state of the industrial automationprocess to determine a first control trajectory of a plurality of inputvariables of the industrial automation process; perform, using thecontrol system, a second search during the predetermined sample periodbased at least in part on the one or more sensor samples indicative ofthe state of the industrial automation process to determine a secondcontrol trajectory of the input variables of the industrial automationprocess; select the first control trajectory for instructing controlactions when the first control trajectory is determined within a desiredtime and is an optimal control trajectory; and select the second controltrajectory for instructing control actions when the first controltrajectory when the first control trajectory is not selected.
 31. Thecontrol system of claim 30, wherein the processing circuitry comprises amulti-core processor, wherein: a first core of the multi-core processoris configured to perform the first search, the second search, or both;and a second core of the multi-core processor is configured to controlperformance of the industrial automation process.
 32. The control systemof claim 30, wherein the processing circuitry comprises a multi-coreprocessor, wherein: a first core of the multi-core processor isconfigured to perform the first search; and a second core of themulti-core processor is configured to perform the second search.
 33. Thecontrol system of claim 30, wherein, when the first control trajectoryis not feasible, the processing circuitry is configured to execute theexecutable code to: determine a first stabilizing control trajectory byprojecting the first control trajectory into a feasible region;determine a second stabilizing control trajectory by shifting andpadding a third control trajectory determined in a previous sampleperiod; and determine a first cost associated with the first stabilizingcontrol trajectory and a second cost associated with the secondstabilizing control trajectory based at least in part on an objectivefunction; select the first stabilizing control trajectory as the secondpoint used to initiate the second search when the first cost is lowerthan the second cost; and select the second stabilizing controltrajectory as the second point used to initiate the second search whenthe second cost is lower than the first cost.
 34. The control system ofclaim 30, wherein the processing circuitry is configured to execute theexecutable code to determine that the first control trajectory isfeasible when each control action of the first control trajectory meetseach constraint on performance of the industrial automation process. 35.A method of controlling an industrial automation process, comprising:performing, using a control system, a first search for an optimal firstcontrol trajectory beginning in an infeasible region based at least inpart on one or more sensor samples indicative of a state of theindustrial automation process; performing, using the control system, asecond search for a second control trajectory beginning in a feasibleregion based at least in part on the one or more sensor samplesindicative of the state of the industrial automation process; selectingthe first control trajectory for instructing control actions when thefirst control trajectory is determined within a desired time; andselecting the second control trajectory for instructing control actionswhen the first control trajectory when the first control trajectory isnot determined within the desired time.
 36. The method of claim 35,wherein the desired time comprises a predetermined sample period for theindustrial automation process.
 37. The method of claim 35, wherein thefirst search comprises an infeasible search beginning at a first pointoutside a feasible region for control trajectories.
 38. The method ofclaim 35, wherein the second search comprises a feasible searchbeginning at a second point inside the feasible region for controltrajectories.
 39. The method of claim 35, wherein the first searchconverges on an optimal control trajectory faster than the secondsearch.
 40. The method of claim 35, wherein the first and secondsearches are run sequentially in the sample period.